rm(list=ls()) 
library(data.table)
library(stargazer)
library(lfe)
library(dplyr)
library(tidyr)
library(haven)
setwd("C:/Users/mg5797/Documents/Pollution/ReplicationPackage")


# import temperature: 
tempdf=read.csv("Derived Data/weatherdf.csv",header=TRUE,stringsAsFactors = FALSE)

tempvarsavg="ay_mean_snow  +aymean_high_num_snow+aymean_severe_num_snow+ay_avg_mean_max + ay_avg_mean_min  + ay_avg_bint80  + ay_avg_bint90  + ay_avg_bint100  + ay_avg_bint20 + ay_avg_bint10+ ay_avg_bint0+ ay_mean_prcp +aymean_high_num_prcp+aymean_severe_num_prcp"

districtcovariates="pct_hisp_2005:factor(year) + pct_black_2005:factor(year) + pct_male_2005:factor(year) + pct_white_2005:factor(year) + enroll_2005:factor(year) + pct_ell_2005:factor(year) + pct_speced_2005:factor(year) + pct_FRPM_2005:factor(year) + missing_pct_ell_2005 + missing_pct_speced_2005 + missing_pct_FRPM_2005 "



acs="pctemploy2005:factor(year)  + pctsinglemother2005:factor(year) + pctmanufac2005:factor(year) + pctutility2005:factor(year) + pctbachelorhigher2005:factor(year)"
# import acs

acsdf=read.csv("Derived Data/acs2009_2018.csv",stringsAsFactors = FALSE, header=TRUE)
acsdf$pctemploy=as.numeric(acsdf$pctemploy); acsdf$pctlaborforce=as.numeric(acsdf$pctlaborforce); acsdf$pctutility=as.numeric(acsdf$pctutility); acsdf$pctmanufac=as.numeric(acsdf$pctmanufac)
acsdf$leaid=acsdf$LEAID; acsdf$LEAID=NULL


demogvar_grade="perfrl + perblk + perhsp + perwht  + malpct_test + pcttested+ totenrl "

demogvar_dist="perell + perfrl +perspeced"
acstime="unempall + single_momall + pctlaborforce + pctutility + pctmanufac"

######################### Main Exhibits 
#40km
distance=40
filetoimport=paste("Derived Data/dfbartik2005_",distance,"kmanalysis_annual.csv",sep="")

dfbartik=fread(filetoimport)

dfbartik=merge(dfbartik,tempdf,by=c("leaid","year"),all.x=TRUE)


dfbartik=dfbartik %>% drop_na(cs_mn_all,avgpm25_9, cs_mn_all_lag_cohort) # same drop conditons as last sample
nrow(dfbartik)


# drop missing
dfbartik=dfbartik %>% drop_na(cs_mn_all,avgpm25_9, cs_mn_all_lag_cohort,pctutility2005,pctemploy2005, pctsinglemother2005,pctmanufac2005, pctutility2005, pctbachelorhigher2005, pct_hisp_2005,pct_male_2005,enroll_2005,pct_ell_2005,pct_speced_2005,pct_FRPM_2005, pct_black_2005, pct_white_2005)

nrow(dfbartik)

# instrument: 
instname=paste("bartikleave_",distance,sep="")
# weights:
weightval=paste("totalprod_fuel",distance,sep="")


##### Table 1 
source("Programs/Results/Subprograms/1_Table1.R")

#####Table 3 
source("Programs/Results/Subprograms/2_Table3.R")

source("Programs/Results/Subprograms/TableA7.R")

#####Figure 2
source("Programs/Results/Subprograms/2_Figure2a.R")

setwd("..")
source("Programs/Results/Subprograms/2_Figure2b.R")

setwd("..")
# Rotemberg numbers 
source("Programs/Results/Subprograms/6a_rotemberg_avg.R")